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Abstract 

In this work, we employ a kinetic theory based approach to predict the hydrodynamic forces on elec- 
tromechanical resonators operating in gaseous media. Using the Boltzmann-BGK equation, we investi- 
gate the influence of the resonator geometry on the fluid resistance in the entire range of nondimensional 
frequency variation < ra < oo; here the fluid relaxation time r = jj,/p is determined by the gas viscos- 
ity /j, and pressure p at thermodynamic equilibrium, and u> is the (angular) oscillation frequency. Our 
results support the experimentally observed transition from viscous to viscoelastic flow in simple gases 
at tu> » 1. They are also in remarkable agreement with the measured geometric effects in resonators in 
a broad linear dimension, frequency, and pressure range. 

1 Introduction 

Electromechanical resonators with linear dimensions in the nanometer to micrometer scales are being devel- 
oped for technological applications and fundamental research. One of the most important attributes of these 
nano/microelectromechanical systems (N/MEMS) resonators is that they have very small intrinsic dissipa- 
tion of energy, quantified by their high quality factors Q = O(10 2 — 10 4 ). N/MEMS resonators are thus 
ultrasensitive to external perturbations enabling important technologies such as atomic force microscopy 
(AFM) [2] and bio-chemical sensing [l(J. 

Some of the most promising applications of N/MEMS, however, require their immersion in fluid media 
(e.g., air mixtures or biological fluids), where fluid-device interaction can significantly degrade the overall 
sensitivity [l|, [21], [^(J. Numerous efforts are currently underway to overcome this difficulty and develop 
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future N/MEMS for promising nanotechnological and biomedical applications. Unquestionably, the flows 
generated by N/MEMS demand a novel understanding of fluid mechanics at increasingly smaller time and 
length scales. Conversely, experimental characterization and numerical analysis of fluid-immersed N/MEMS 
resonators provide an invaluable opportunity to advance knowledge in new areas of fluid mechanics such as 
high-frequency nanofluidics. Recent work on high-frequency oscillating flows 
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reports a viscoelastic 



transition in simple gases at sufficiently large values of the nondimensional frequency tlo. Here, r = /i/p is 
the relaxation time in terms of the pressure p and viscosity /i of the gas at equilibrium; ui is the oscillation 
frequency. Such remarkable phenomenon is beyond the reach of classical (Newtonian) fluid mechanics, which 
is only valid for tlo <C 1. The viscoelastic transition is accompanied by a substantial attenuation of the energy 
dissipated by the fluid and a subsequent improvement in the performance of the fluid-immersed device. 

In the present study, we solve the Boltzmann-BGK equation of kinetic theory via appropriate numerical 
procedures 8|, |23j in order to predict fluidic effects (e.g., damping force and energy dissipation) for specific 
geometries. The validity of the kinetic methods is not constrained to Newtonian flow; thus, our models yield 
excellent agreement with experimental measurements on different resonators over a wide range of frequency 
and pressure variation. 



2 Electromechanical resonators 

Illustrated in figures HJa)-[lJb) are the first class of studied resonators in the form of cantilever and doubly- 
clamped beams (L z 3> L x ~ L y ). Harmonic motion in the beams can be induced through the application of 
periodic electrostatic, photothcrmal, or incrtial forces. The beams are driven around their fundamental and 
first harmonic out-of-plane flexural resonances while optical techniques are used to determine the resonant 
response [16| . The doubly-clamped beams are suspended above a stationary substrate at a mean height 
A ~ 400 nm; thus, the presence of the substrate has no significant effect on the fluidic damping (e.g. via 
squeeze-film damping) [15| . We also study a macroscopic quartz crystal resonator [11| . The studied resonator 
[see figure (He)] is in the form of a thin crystal disk (L x = L z ^> L y ) connected to electrodes so that its 
resonances in thickness-shear modes can be electrically excited and detected by piezoelectric effects. In all 
the measurements, the resonance amplitudes of the beams and crystal resonators are kept extremely small. 

Specific dimensions and (vacuum) characteristics, such as resonance frequency u> , quality factor Q Q , and 
surface to modal mass ratio S/ m , of four studied resonators are listed in table O The size of the devices vary 
from sub-micron to millimeters while their resonance frequencies are in the range of kilohertz to megahertz. 
Experimental analysis of the four resonators in table [2] is performed with the devices operating in purified 
nitrogen at room temperature T ~ 300 K; the pressure is gradually varied from low- vacuum to atmospheric 
pressure, 0.1 < p < 1000 Torr. As the pressure is varied and given that r = fi/p, the resulting flows cover a 
wide range of dimensionless frequency variation 0.001 < tlo < 10 [lsj . 
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Class 

(*: first harmonic) 






Dimensions 




U> /2lT 

(MHz) 


Qo 


S/m 
(m 2 /kg) 


AR 


Cantilever Beam 


L x 


= 2.0 


Ly — — 53 L z = 


460 


0.078 


8320 


688 


26 


Cantilever Beam* 


L x 


= 3.6 


Ly = 36 L z = 


125 


1.97 


3520 


405 


10 


Doubly-Clamped Beam 


L x 


= 0.2 


L y = 0.23 L z = 


9.6 


24.2 


415 


8380 


1 


Quartz Crystal Disk 


L x 


=L Z = 


= 10000 Ly ~ 


160 


32.7 


40755 


6.29 






Table 1: Electromechanical Resonators 
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Resonators immersed in a fluid. The dynamics of an electromechanical resonator immersed in a fluid can 
be studied by means of a one-dimensional harmonic oscillator approximation Q 

m {x tt +j x t +Lu 2 x) = F(t) + F f (t) (1) 

where m a is the effective mass corresponding to the vibrational mode HQ , 7„ is the (structural) damping 
coefficient and oj is the resonance frequency of the device in vacuum. The load on the oscillator is produced 
by the driving force F(t) — Re{F(uj)e~ luJt } along with a fluid resistance Ff(t) = — m G (7/ xt + (3f Xu) 



which has both dissipative and inertial components [lTj. The oscillation amplitude has the general form 
x(t) = Ke{x(u>)e~ l ^ t+ ^} and from equation fl} the frequency response is 

X (u)) = — — 1 — ■ -r, (2) 

TO [UJ Z -UJ\1 -lw( 7 o+7/)J 

with Fq = F(u})e 1 ^ being the effective force amplitude. Equation @ includes fluidic effects through the 
fluidic inertia, or fluid-added mass, j3f and fluidic damping 7/. Experimental values of the added mass and 
fluidic dissipation are respectively obtained from the frequency shift and the broadening of the Lorentzian 
frequency response ([2]) [15| . In this study, the fluid-added mass is very small /3j < 1 and so is the measured 
shift in the resonance frequency: Alo /lu w 0f/2. The device quality factor in the fluid thus becomes 

Q~UJo/{lo+lf)- 

3 Hydrodynamics of high-frequency flows 

Similar to previous work [l , [21], , the present analysis is valid when gradients along the spanwise direction 
of the oscillating body are negligible so that the flow is considered two-dimensional. This assumption holds for 
a slender beam (L z 3> L x ~ L y ) with uniform rectangular cross-section, or a thin disk (L x = L z = tf 3> L y ), 
where the aspect ratio AR = L y /L x (see table [2]) becomes the geometric parameter characterizing the 
dominant chordwise flow. Other dynamically relevant flow parameters are the Mach number M — U /c s 
and the Reynolds number Re = pU Ly/iJL determined by the fluid velocity amplitude U = \u)x{ui)\, the 
speed of sound c s ~ 298 m s _1 , and molecular viscosity /i ~ 1.78 x 10 -5 kg m" 1 s _1 of nitrogen at room 
temperature. Given the tiny oscillation amplitude of the resonators, the resulting Mach and Reynolds number 
are extremely small, M ~ Re < 0.001 in all operation conditions. Consequently, the generated flow can be 
assumed laminar, nearly incompressible, and isothermal. 

The common approach hitherto encountered in the literature [see[l|, [21} [2^] for the theoretical and/or 
numerical determination of fluidic effects on resonators is based on the Navier-Stokes equations for Newtonian 
fluid flow. However, a fundamental assumption for the applicability of classical Navier-Stokes equations is 
that hydrodynamic scales are much larger than their kinetic counterparts, i.e. tlj <C 1. Previous work 
[3, Q, S, U I on high-frequency oscillating flows has established that the non-dimensional frequency tlo 



4 



determines qualitatively different behavior: purely viscous (Newtonian) for tlo = 0, viscoelastic (transitional) 
for < tlo < oo, and purely elastic (free- molecular) for tlo — > oo. While viscous flows are accurately described 
by Newtonian hydrodynamic equations, an extended hydrodynamic description accounting for kinetic (non- 
equilibrium) phenomena is required at sufficiently large tlo where non-Newtonian behavior is observed [see 
26, \ 



3.1 Newtonian hydrodynamics 

In the Newtonian regime tlo -C 1, the flow around a body oscillating with very small amplitude |x(a;)| <C L 

m 



is governed by the linearized Navier-Stokes (NS) equations for incompressible flow 

V-u = 0, — = vV 2 u~ -Vp. (3) 
at p 

After adoption of the standard no-slip boundary conditions, these equations provide well-known analytical 
solutions for simple geometries [17I Yl\ . Through these solutions, one can determine fluidic forces over a 
body as 

F f (t) = pBLo 2 Ke{T{Lo)x{Lo)e- l ^ t+ ^}, 

where p is the fluid density and B is the body volume; T = m (j3f+i^ff /uj)/(pB) is the so-called hydrodynamic 
function. One of the simplest solutions of the unsteady NS equations is obtained for an infinite plate 
oscillating in a fluid (L x = L z = oo,L y = 0), known as the Stokes' second problem 26, For a slender 
body having a thin cross section with small but finite width (0 < AR <C 1) the Newtonian hydrodynamic 
function can be approximated by the solution of the Stokes' second problem: 

r (ylH<1) ( w ) = (i + <) - (4) 

Here, S is the surface area in contact with the fluid and v = fi/p. When the cross section is not small 
{AR > 1), the common approach has been to study the flow generated by simpler geometries such as a 
cylinder with a radius equal to half the nominal length scale, r = L y /2. For the case of an infinite cylinder 



with its axis normal to the ^-direction, the hydrodynamic function is 



1Q,Q 



r ! ,\ 1 j Ai K^-iy/iLor 2 /^) 

I cyl(U) = H : -, (5) 

^JiLor 2 /v ¥±o(—i^/iLor 2 /v) 

where Ko and K\ are Bessel's functions of the third kind. The asymptotic behavior of the Newtonian 
hydrodynamic function of a cylinder {5} and that of a rectangular beam are similar in both limits lot 2 jv — > 
and lot 2 jv — > oo [24]. For this reason, and only within the Newtonian regime (tlo <C 1), the solution 
for an oscillating cylinder has been employed with some degree of accuracy in estimating fluidic effects 
on rectangular beams with cross-sections of moderate to large aspect ratio [l|, [21I 20 1. For very large 



aspect ratios (L x — > 0), 21( formulated an empirical correction to equation (JSJ such that the hydrodynamic 
function becomes — T cy i£l(uj L 2 / v) . Nevertheless, the correction of [21] remains essentially unity 

|fi — 1| = O(10 _1 ) within all regimes studied in this work. 



■5 



3.2 Beyond Newtonian hydrodynamics 

When flow time scales T = l/u> are of the same order as the relaxation time r, kinetic effects become 
significant and Newtonian hydrodynamic equations ([3]) breakdown. The primary issue encountered beyond 
Newtonian regimes (rw > 1) is the lack of a robust hydrodynamic equation governing the flow. To obtain 
a hydrodynamic description valid for arbitrary nondimensional frequencies (0 < tlj = t/T < oo), one must 
resort to kinetic theory representations of the flow. Unfortunately, the problem of deriving (closed-form) 
hydrodynamic equations via kinetic theory, albeit largely studied [see Q Q, S, 3] , remains essentially open 
for arbitrary flow regimes. 

Another critical point arising when kinetic effects are no longer negligible is that of proper boundary 
conditions for the hydrodynamic equations at the solid-fluid interface. Hydrodynamic boundary conditions 
are determined by a rather complex fluid-solid interaction in the vicinity of a solid surface. For flows at 
finite Knudsen number, the slip boundary condition has been extensively adopted 0, Q, According to 
Maxwell's picture of slip of a gas over a solid surface, a finite mean-free-path A ~ tc s leads to an effective 
slip velocity 

[u(x B , t) - XJ W ] ■ t = — — A [V(u • t) • n + V(u • n) • t] (6) 

a v 

to be employed as boundary condition at the coarse-grained (hydrodynamic) level. Here t and n are the 
unit tangent and normal vectors to a wall located at x„ and moving with velocity \J W ; u is the fluid 
velocity. Meanwhile, a v is the tangential momentum accommodation coefficient of the solid surface. A unit 
accommodation coefficient (a v — 1) represents a situation where all fluid particles are diffusively scattered 
after collision with the wall; the opposite limit (a v = 0) corresponds to the case where all such collisions are 
specular. First-order Maxwell slip models ^ are accurate for steady flow at small to moderate Knudsen 
numbers (Kn < 1) [2, [ill while high-order versions have been proposed for unsteady shear flow QjQj- 
The effective slip in oscillating shear flows has recently been investigated via kinetic methods such as lattice 
Boltzmann-BGK (LBGK) and direct simulation Monte Carlo (DSMC). Expected functional shapes have been 
obtained for the slip as a function of the Knudsen number Kn — X/C in steady shear flows (£ ~ | Vut|/|ut|) 
or the nondimensional frequency tlo = t/T in oscillating shear flows 



4 Kinetic model of hydrodynamics 

At room temperature and under ordinary pressure conditions (ranging from low-vacuum to atmospheric 
pressure), simple gases are composed of a large number of electrically neutral molecules, each with an 
effective diameter that is negligible compared to the average intermolecular distance. Under such conditions 
the Boltzmann-BGK equation (BE-BGK) is an accurate kinetic model of the flow [4|,[6|]. For monatomic gases 
in the absence of external force fields F = 0, the BE-BGK for the evolution of the Boltzmann distribution 



G 



/(x, v, t) in phase space (x, v) reads: 



df f - f eg 

-L+V.Vf = -± J —. 

at t 



(7) 



Without loss of generality, we define 9 = IcbT '/m gas = c 2 and adopt a molecular mass m gas = 1; the 



equilibrium distribution can then be expressed as 



r 9 (x,v,i) 



(2ir9) : 



■ exp 



(v-u) 2 



29 



(8) 



where D is the velocity space dimension (v = Ufce^; k = 1, -D). Hydrodynamic quantities, like mass density 
p, fluid velocity u, and energy are obtained as moments of the distribution function: 



J /(x,v,t)dv = p(x,t), 
J f{x,v,t)vd\ = pu(x,t), 
J /(x,v,<)v 2 dv = pD0(x,t) +/m 2 (x,<). 



(9) 



Kinetic boundary conditions. For bounded flows, particular solutions of equation ([7]) will require proper 
boundary conditions. Within the framework of classical kinetic theory we consider the gas as bounded by a 
perfectly elastic and isothermal surface (9 W = 9) located at x m while moving with velocity \J W . Under this 
depiction, particles impinging on a solid surface with velocity v' acquire a post-collision velocity v defined 
by the scattering kernel £>(v' — » v). General boundary conditions will then read [4]: 

\{v-V w )-n\ f(x w ,v,t)= f\(v'-V w )-n\ B(v' ^ v)/(x w , v', f)dv'; (10) 

(v'-U ro )-n<0 

for (v — TJ W ) • n > 0. In this work we implement and assess two different kinetic boundary conditions at the 
fluid-solid interface; diffuse scattering (DS) 

(v-U lu ) 2 " 



/K,v,t) 



(2ir9y 



■ exp 



29 



(v - U w ) • n > 0, 



and bounce-back (BB) 



/(x™,v,f) =/(x w ,-v + 2U w ,t); (v-U u ,)-n>0, 



(11) 



(12) 



where impinging particles [(v — U m ) • n < 0] are reflected back with the same relative speed and angle of 
incidence. While a diffuse-scattering (DS) kernel yields hydrodynamic slip in agreement with the Maxwell 
slip model ^ for a fully accommodating surface (a v = 1), the bounce-back model renders no-slip at the wall 
for all flow regimes (i.e. the surface has negative accommodation coefficient a v = —2). The validity range of 
each model for the solid-fluid interaction will be assessed when comparing LBGK simulation using both DS 
and BB schemes against experimental results. 
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4.1 Free-molecule hydrodynamics 

The nondimensional frequency tlu = 2tt\/Cfm is proportional to the ratio of the (equilibrium) mean free 
path A = rc s to the average distance Cfm = c s 2tt/lu traveled by a particle during one oscillation period. 
Hence, in the limit tlu — ► oo, the kinetic dynamics becomes practically collisionlcss and free-molecular flow 
approximations are applicable. Let us now analyze a L x x L y rectangular section in free molecular flow 
moving in the x-direction at velocity u = ui such that M = U w /V6 C 1. In the case of diffusive wall 
scattering (fTTjl the net x-force per unit length is Q] 



Ff s (t) = 




For the bounce-back model (|12p . which renders no hydrodynamic slip, one has Ff B = 2Ff s . Similar 
expressions can be obtained for other models of the gas-surface interaction, e.g. Fj P {t) = — (iJ~^Ly \ p\fd u 
for specular reflection models rendering no shear stress. Clearly, in free molecular flow (tlu — > oo), there is 
no fluidic inertia (/?/—> 0) and hydrodynamic forces only have a dissipative component (Ff = fd it)- 



5 Lattice Boltzmann BGK simulation 



23 



22] and 



The method in this work falls in the class of high-order LBGK models originally formulated by 
investigated by Q for isothermal unidirectional flow in non-Newtonian regimes. Discretization of velocity 
space in a finite number of lattice velocities {v^; i — 1,Q} allows one to reduce the problem of solving 
equation (O to that of solving a set of lattice Boltzmann-BGK (LBGK) equations 

^+v i .V/ i = -^^ (14) 

Ot T 

jf =Ml + 5(v,'u) + ^(vru) 2 -y] i = l,Q (15) 

governing the evolution of each lattice population fi(x,t) = i«i/(x, v;,t) in configuration space. The set of 
LBGK equations (|14M15[) is formally derived by projecting the continuum Boltzmann-BGK equation ([7][5]) 
onto the Hilbert space H 2 spanned by the orthonormal basis of Hermite polynomials up to second order. 
A Gauss- Hermite quadrature formula determines the lattice velocities v i5 i.e. integration points, and their 
associated weights Wi. A quadrature rule with algebraic degree of precision d > 4 permits the exact numerical 
integration of the three leading moments (|10[) of the continuum distribution /: 

Q Q Q 

p(x,t) = 5^/i(x,i), pu(x,t) = 5^/ i (x,t)v i , P (u 2 + D9){ x ,t) =£/ i (x,t)v?. (16) 

i—l i—1 i—1 



The particular lattice employed for the present simulations is the D2Q37 [22j (i.e. velocity space dimension 
D = 2, number of lattice velocities Q = 37); the weights and velocity abscissae of the lattice are included 
in the appendix. After lattice discretization of configuration space (Ax, = v^A<) numerical procedures to 
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solve the LBGK equations (fH| advance in two steps: advection and collision. During the advection step the 
streaming of lattice populations {/j; i = 1, Q} is performed along the lattice directions: 



/^(x,t) = / i (x-v i At,t). 



(17) 



At the end of the advection step, mass, momentum, and energy (|16[) are computed using f? dv allowing the 
explicit evaluation of equilibrium distributions (1151) . Then the BGK ansatz is applied at the collision step: 



/ i (x + v i ,t + At) = /? fl + 
The projected, or dealiased, nonequilibrium component 



T 



Si 



27] 



sr = ^{v aiV0% - es aP ) Y^uf v - s; q )v aj v 0j a ,p = i,D % = 1, q 



(18) 



(19) 



3=1 



ensures that the post-collision distribution /j is strictly expressed by a linear combination of up to second- 
order Hermite polynomials, i.e fi £ H 2 . The regularization procedure [27} specified in equations (|19[) - (|18[) 
eliminates significant numerical errors due to lattice orientation anisotropy and aliasing effects that are 
reported for the standard LBGK algorithm in strong non-equilibrium conditions |s , 27 1 . 



6 Results 

The quantitative analysis of fluidic damping on devices with different sizes and structural features is per- 
formed by employing the nondimensional quantity: 

j f (roj,AR)= 71 . (20) 

Here, S/m is the ratio of wet area to modal mass reported in table [2 After normalization by the fluid 
mass density p and speed of sound c s = >/6, the dimensionless damping 7* (I20p is solely determined by 
the nondimensional frequency tlo and hydrodynamic shape characterized by the aspect ratio AR — L y /L x . 
Based on reported data in table [21 we adopt AR ~ for the quartz resonator and AR ~1, 10 & 26 for the 
beams. Experimentally measured values of the total damping are presented as open symbols in figure [2] 
These are determined from measurements of the resonant response © as a function of pressure as discussed 
above. As observed in figure [2l experimental values of the density-normalized damping 7^ (|20p exhibit an 
increasing dispersion in the upper limit of dimensionless frequency variation tlo — > 00 where 7/ — * 0. This 
is due to the fact that the measured fluidic effect becomes smaller and smaller as the pressure is lowered 
and the measured signal is dominated by the finite error ~ 5% present in all the experimental data. This is 



discussed in more detail in 
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Lattice Boltzmann-BGK (LBGK) simulation. Simulations are performed with the D2Q37-H2 model 
described in Sec. \5\ It has been demonstrated that high-order LBGK models such as D2Q37-H2 with 
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a regularization procedure yield excellent agreement with extended hydrodynamic descriptions derived for 
isothermal and unidirectional flows in both Newtonian and non-Newtonian regimes [si, [3]. In all simulations, 
t = fi/p = v/9 is readily determined from the gas properties at thermodynamic equilibrium, while the 
employed Mach number is very low M = U o /V0 < 0.01. Simulation results on different resonators obtained 
by the D2Q37-H2 model for tuj— 0.001, 0.01, 0.1, 0.5, 1, and 10 are compared against experimental data 
and available analytical expressions for Newtonian and free-molecular flow in figure [5] The employed LBGK 
models yield a remarkable agreement with experimental measurements over a wide range of nondimensional 
frequency 0.001 < tuj < 10 for different device geometries and dimensions (table [2]). On the other hand, 
Newtonian fluid approximations such as equation (U) for the quartz disk (AR ~ 0) and ([5]) for the beams 



(AR > 1) [l,|2l|,|20| give acceptable agreement only within the low frequency limit tuj < 0.1. It is important 



to remark that all Newtonian predictions largely overestimate the fluidic dissipation in the frequency range 
tuj > 0.1. In the high-frequency limit tuj 1, only kinetic approaches such as LBGK simulation and the 
free- molecular flow solution given by equation (|13[) are in good agreement with experiment. 

Fluid-solid interaction and boundary schemes. Owing to the small amplitude of oscillation (|a;(ii;)| <C L x ), 
the solid boundary can be assumed to remain static. Only the boundary velocity U w = U a sin(ujt)i varies in 
time with velocity amplitude U a < ujAx such that the displacement amplitude is smaller than the distance 
between neighboring lattice nodes. Two boundary schemes are employed in LBGK simulation when modeling 
the moving walls: the diffuse-scattering scheme (DS) and the bounce-back scheme (BB) that were explained 
in Sec. |U While the BB scheme (fT2j) renders no-slip as hydrodynamic boundary condition, the DS model 
(fTTjl produces an effective slip velocity between the body surface and the fluid immediately adjacent to 
it. Experimental results in figure [2] show that numerical schemes rendering no-slip at the wall overpredict 
the fluidic dissipation when tuj > 0.1. On the other hand, the diffuse-scattering scheme which models the 
resonator surface as fully accommodating (a v = 1) yields a close agreement with experimental data in the 
entire studied range 0.001 < tuj < 10. 

Quality factor and fluidic effects. Since fluidic inertia is very small ([3f <C 1) for the studied flow conditions 
the quality factor of the fluid-immersed device is 

Q = Qo j= • (21) 

Evidently, dissipative effects quantified by the structural dissipation 7 D = uj /Q and fluidic damping 7^ = 
P\j\ ^~7/ determine the device performance. The quality factors reported for the four studied devices 
are compared against numerical predictions in figure [3] We observe that lattice Boltzmann-BGK simulation 
(model D2Q37-H2 with DS wall treatment) is in close agreement with experimental data obtained for different 
cross-sections (0 < AR < 26) in a wide range of operation conditions (0.1 < p < 1000 Torr, 0.078 < uj < 
24 MHz). 
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Figure 2: Nondimensional fluid damping: ^^(tuj,AR) vs. to. (a) uj /2-k = 32.7 MHz, 
AR = 0. (b) w /27r = 24.2 MHz, AR=1. (c) w /2tt = 1.97 MHz, AR=10. (d) cj o /2vr = 0.078 
MHz, Ai? = 26. Open symbols (red): experimental data. Closed dots (black): LBGK 
simulation rw=0.001, 0.01, 0.1, 0.5, 1, and 10; solid line (black): D2Q37-H2 DS; dotted 
line (black): D2Q37-H2 BB. Dashed line (blue): Newtonian fluid approximations. Dashed- 
dotted line (green): free molecular flow (DS). 
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Figure 3: Quality Factor: 



10 



10 



10 



7o 



7o+7/ 



vs. p [Torr]. Open symbols (red): experimental 



data; (O) oj /2ir = 32.7 MHz, AR = 0; (O) w /2tt = 24.2 MHz, AR=1; (A) w /2tt = 1.97 
MHz, AR=10; (□) w /2tt = 0.078 MHz, AR = 26. Closed dots (black): LBGK simulation 
(D2Q37-H2 DS) rw=0.001, 0.01, 0.1, 0.5, 1, and 10. 

7 Conclusions 



Experimental measurements of fluidic effects on diverse electromechanical resonators have been compared 
against available analytical approaches, i.e. the Newtonian flow and free-molecular flow approximations, 
as well as against kinetic-based simulations presented in Sec. [5] Clearly, Newtonian and free-molecular 
flow models fail to describe transitional flow in the region 0.1 < tlo < 10. On the other hand, lattice 
Boltzmann-BGK (LBGK) simulation accounting for specific geometrical features accurately represents the 
fluidic damping 7/ in all studied regimes: Newtonian (™ < 0.1), transitional (0.1 < tlo < 10), and free- 
molecular flow (tlo > 10). 

Viscoelastic dynamics in high-frequency flows. The invalidity of Newtonian approaches as tlo — > 00 is 
not only due to surface effects, which might be absorbed by prope r hy drodynamic boundary conditions, but 



also due to the qualitatively different fluid dynamics in the bulk [26, 8]. Kinetic-based (LBGK) simulation 
precisely reproduces the experimentally observed (figure [2]) saturation of density-normalized dissipation 
(7 j — > const) in the high-frequency limit tlo — > 00. This remarkable phenomenon involves a gradual 
transition from viscous to viscoelastic to purely elastic flow of a simple monatomic gas that has been reported 



by previous theoretical 



29 . l8( and experimental [15:] studies. The viscoelastic response of simple gases in 



the high-frequency limit is a well-known phenomenon within the realm of transport theory and statistical 
12j | ; diffusion processes in nonequilibrium systems can only be established after a finite time ~ r 



physics 



of the order of the relaxation time. In the short-time limit t < t, where diffusion effects are still weak and 
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12| , one observes an ensuing decay 



transport coefficients such as shear viscosity become frequency dependent 
in the dissipation of fluid momentum and energy. 

Near-wall phenomena. Kinetic effects can never be neglected within the so-called Knudsen layer, e.g. at 
distances from the solid boundary that are smaller than one mean free path. Thus, kinetic parameters, such 
as the relaxation time and mean free path, are not easily determined in the near-wall region where gas-solid 
interaction is significant. In consequence, within a mean free path from the wall, Newtonian fluid models 
for the stress break down and NS equations must actually be applied immediately outside the Knudsen 
layer. The concept of having effective slip as proper boundary condition for hydrodynamic (coarse-grained) 
equations must be understood within this context. Our kinetic model of the flow based on the Boltzmann- 
BGK is rather simple and does not accurately represent the Knudsen layer; it relies on a single constant 
relaxation time r = fi/p and kinetic boundary conditions (I10|) determined by surface scattering kernels 
B(y' — > v) for a perfectly elastic and isothermal surface. Nevertheless, the kinetic model in this work 
accurately predicts hydrodynamic effects such as fluid resistance and mean energy dissipation via adoption 
of a Maxwell scattering kernel with surface accommodation a v — 1 (i.e. the DS scheme explained in Sec. [4}. 
The net effect of the studied gas-surface interactions can be assessed by comparing results (figure [5]) from 
DS (slip) and BB (no-slip) schemes in the range 0.001 < tuj < 10 for different cross-sections < AR < 26. 
The reduction in fluid damping solely due to effective slip is found in the interval — 0.4. 

Resonator performance beyond Newtonian regime. The decay of the energy dissipated by the fluid as 
tuj — > oo has beneficial effects on the resonator performance. Under relevant experimental conditions, 
fluidic damping largely dominates over structural dissipation, 7/ 3> 70, and thus, Q « u/jf. In such 
conditions the quality factor will actually increase linearly with the operation frequency, Q cx u> for tuj > 1, 
instead of the square root dependence, Q cx y/uj, observed for Newtonian flow tuj <C 1. Therefore, it is 
advantageous to operate the resonator at a frequency tuj > 1 well within the viscoelastic regime. This 
could be accomplished either by increasing the resonance frequency or by decreasing the effective relaxation 
time of the fluid through less trivial mechanisms, e.g. through polymer addition or foams for water. Other 
strategies that can potentially improve the device performance include modifying the cross-sectional shape 
and surface properties. As observed in figure[2]for the entire range 0.001 < tuj < 10, bodies with small aspect 
ratios AR = L y /L x <C 1 generate less fluidic damping at the same operation frequency and surrounding gas 
conditions. On the other hand, the employment of superhydrophobic coatings for resonators in water can 
further increase the effective hydrodynamic slip with a subsequent reduction of the resistance forces. 

Lattice Boltzmann-BGK simulation for N /MEMS Hydrodynamic. The qualitative and quantitative agree- 
ment between our LBGK simulations and experimental data over a wide range of pressure 0.1 < p < 1000 
Torr and frequency variation 0.001 < tuj < 10 constitutes a remarkable achievement for the kinetic method- 
ology applied in this work. The precise determination of fluid forces and, thus, device quality factors for 
diverse geometrical configurations in widely different operation regimes is of fundamental importance in 
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states 




(0,0) 


1 


0.233150669132352000 


(1,0) 


4 


0.107306091542219000 


(±1,±1) 


4 


0.057667859888794800 


(2,0) 


4 


0.014208216158450700 


(±2, ±2) 


4 


0.001011937592673570 


(3,0) 


4 


0.000245301027757717 


(±1,±2) 


8 


0.005353049000513770 


(±1,±3) 


8 


0.000283414252994198 



"Lattice Speed: c = 1.19697977039307 \fd 
Table 2: Model Parameters D2Q37 
advancing the design of future N/MEMS devices. 
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A Appendix: The LBGK model D2Q37 

Velocity abscissae and weights of the D2Q37 lattice model 



23, 



22J is presented in table [2] 
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